############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(sf)
base::library(sp)
base::library(colorspace)

drug_data <- read.csv(here("Data", 
                           "NCHS_-_Drug_Poisoning_Mortality_by_County__United_States.csv"))
glimpse(drug_data)
table(drug_data$Year)

deaths_2004 <- drug_data %>% filter(Year == 2004)

deaths_2020 <- drug_data %>% filter(Year == 2020)


county_sp <- st_read(here("Data", "shapefile",
                          "cb_2018_us_county_500k.shp"))
glimpse(county_sp)
class(county_sp)
table(county_sp$STATEFP)

county_sp <- county_sp %>%
  mutate(GEOID = as.numeric(GEOID)) %>%
  mutate(STATEFP = as.numeric(STATEFP)) %>%
  select(FIPS = GEOID, FIPS.State = STATEFP, ALAND, AWATER, geometry)

county_sp <- county_sp %>% #filter(!(state_code == 2)) %>%
  filter(FIPS.State < 60)
table(county_sp$FIPS.State)

data_2004 <- merge(county_sp, deaths_2004, by = c("FIPS", "FIPS.State"))
glimpse(data_2004)

data_2020 <- merge(county_sp, deaths_2020, by = c("FIPS", "FIPS.State"))
glimpse(data_2020)



###### 2004 Map
data_proj04 <- data_2004 %>% 
  filter(!(State == "Alaska")) %>% filter(!(State == "Hawaii"))

data_proj04 <- st_transform(data_proj04, crs = 2163)
glimpse(data_proj04)

summary(data_proj04$Model.based.Death.Rate)
summary(data_2020$Model.based.Death.Rate)

plot04 <- ggplot(data_proj04) +
  geom_sf(aes(fill = Model.based.Death.Rate), colour = NA) +
  ggtitle("A. Death Rates 2004") +
  scale_fill_continuous_sequential(palette = "Reds",#, #trans = "reverse",
                                   name = NULL, limits = c(0, 143),
                                   na.value = "grey50") +
  coord_sf(crs = st_crs(data_proj04), datum = NA) + 
  theme_minimal() + ylab(NULL) + xlab(NULL) + theme(legend.position = "right") +
  theme(title = element_text(size = 13))

plot04




###### 2020 Map
summary(data_2020$Model.based.Death.Rate)

data_proj20 <- data_2020 %>% 
  filter(!(State == "Alaska")) %>% filter(!(State == "Hawaii"))

data_proj20 <- st_transform(data_proj20, crs = 2163)

plot20 <- ggplot(data_proj20) +
  geom_sf(aes(fill = Model.based.Death.Rate), colour = NA) +
  ggtitle("B. Death Rates in 2020") +
  scale_fill_continuous_sequential(palette = "Reds",#, #trans = "reverse",
                                   name = NULL, limits = c(0, 143),
                                   na.value = "grey50") +
  coord_sf(crs = st_crs(data_proj20), datum = NA) + 
  theme_minimal() + ylab(NULL) + xlab(NULL) + theme(legend.position = "right") +
  theme(title = element_text(size = 13))

plot20

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

multiplot(plot04, plot20, cols = 1)





